We  describe  the  construction  of  a  collection  of  quadrature  formulae  suitable  for  the  efficient  dis¬ 
cretization  of  certain  boundary  integral  equations  on  a  very  general  class  of  two-dimensional  do¬ 
mains  with  corner  points.  The  resulting  quadrature  rules  allow  for  the  rapid  high  accuracy  solution 
of  Laplace’s  equation  and  the  Helmholtz  equation  on  such  domains.  Our  approach  can  be  adapted 
to  many  other  boundary  value  problems  as  well  as  to  the  case  of  surfaces  with  singularities  in  three 
dimensions.  The  performance  of  the  quadrature  rules  is  illustrated  with  several  numerical  examples. 
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1.  Introduction 


In  [4] ,  a  procedure  for  the  construction  of  an  orthonormal  basis  spanning  the  space  of  restrictions  of 
functions  satisfying  a  Laplace  boundary  integral  equation  over  a  contour  T  to  a  small  curve  segment 
r0  C  T  is  introduced.  Using  the  machinery  of  [3],  such  a  basis  can  be  exploited  to  form  a  “purpose- 
made”  quadrature  formulae  for  the  efficient  discretization  of  the  boundary  integral  equation  over  the 
curve  segment  IV  As  the  size  of  the  resulting  quadrature  rules  depends  only  weakly  on  the  geometry 
of  the  curve  segment  IV  boundary  integral  equations  over  curves  with  complicated  geometry  (e.g., 
corners)  can  be  discretized  quite  efficiently.  Moreover,  by  repeating  the  basis  construction  procedure 
for  multiple  curve  segments,  a  collection  of  quadratures  suitable  for  the  efficient  discretization  of  a 
particular  boundary  integral  equation  over  an  entire  class  of  contours  with  complicated  geometry  can 
be  precomputed.  We  refer  to  such  a  collection  as  a  set  of  “universal  quadratures.”  The  construction 
of  universal  quadratures  for  boundary  integral  equations  arising  from  Neumann  and  Dirichlet  Laplace 
boundary  value  problems  on  polygonal  domains  is  described  in  [4].  The  resulting  quadrature  rules 
achieve  double  precision  accuracy  with  approximately  30  discretization  nodes  required  per  corner. 

In  the  present  paper,  we  introduce  a  modified  version  of  the  procedure  of  [4]  applicable  to  boundary 
value  problems  for  both  Laplace’s  equation  and  the  Helmholtz  equation  with  small  wavenumbers.  We 
also  describe  the  construction  of  collections  of  universal  quadratures  for  the  corresponding  boundary 
integral  equations  on  a  very  general  class  of  domains  with  corners.  The  resulting  quadratures  allow  for 
the  rapid  high  accuracy  solution  of  Laplace  and  Helmholtz  boundary  value  problems  on  such  domains. 

A  large  number  of  schemes  for  the  solution  of  boundary  integral  equations  on  domains  with  corners 
have  been  proposed,  but  unlike  the  work  presented  here,  these  techniques  all  suffer  from  one  or  more 
of  the  following  serious  disadvantages: 

1.  They  require  dense  meshes  of  discretization  nodes  near  corner  points. 

It  is  possible  to  discretize  boundary  integral  equations  near  a  singular  point  of  the  boundary  by 
simply  using  a  very  dense  quadrature.  This  practice,  however,  results  in  very  large  systems  of  linear 
equations  in  the  case  of  domains  with  many  corner  points. 


2.  They  involve  classical  quadrature  techniques  like  subtraction  of  singularities,  factorization,  or  sub¬ 
stitution. 


These  techniques  are  simply  too  cumbersome  to  apply  to  large-scale  domains.  The  behavior  of  a 
solution  to  a  Helmholtz  boundary  integral  equation  with  wavenumber  k  near  a  corner  point,  for 
instance,  can  be  described  as 


a{x) 


r]{x)  +  a(x,  k)  for  x  <  0 

r](x)  +  /3(x,  k)  for  x  >  0  ’ 


where  ry  is  a  singular  function  determined  by  the  angle  of  the  corner  and  a  and  (3  are  oscillatory 
functions  whose  behavior  is  dictated  not  only  by  the  angle  of  the  corner  point  but  also  the  geometry 
of  the  boundary  curve  on  either  side  of  the  corner  point  and  the  wavenumber  of  the  Helmholtz 
equation.  It  might  be  that  any  particular  case  can  be  handled  using  classical  quadrature  techniques, 
but  the  number  of  modifications  required  in  order  to  apply  these  techniques  to  the  domain  shown  in 
Figure  5c  of  Section  7,  for  instance,  is  already  daunting  and  it  is  by  no  means  a  truly  “large-scale” 
example. 

Perhaps  even  more  problematic  is  the  lack  of  straightforward  generalizations  of  these  approaches 
to  surfaces  with  singularities  in  M3. 


3.  They  require  the  modification  of  the  boundary  integral  equation. 

For  instance,  many  approaches  for  the  solution  of  Laplace’s  equation  on  domains  with  corners  take 
advantage  of  the  fact  that  the  second  kind  boundary  integral  equation  over  the  contour  T 

1  f  d 

±2a(x')  +  Jr  5j7loglx  -  y I  a(x)dS(y)  =  u(x)> 
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where  d/dvy  denotes  differentiation  in  the  variable  y  with  respect  to  the  outward  normal  vector  of 
T,  is  related  to  the  Cauchy  integral 


a(z) 


dz. 


(1.1) 


Jr  z-w 

Operations  involving  the  equation  (1.1)  can  often  be  reduced  to  operations  on  the  more  tractable 
integral 

f  cr(z)  —  cr(w) 

Jr  z  —  w 


-dz, 


which  is  nonsingular  when  a  is  sufficiently  smooth.  Unfortunately,  approaches  of  this  type  do 
not  apply  to  many  boundary  integral  equations  of  interest,  and  even  when  they  do,  they  require 
adaptation  to  each  new  case.  Moreover,  this  class  of  techniques  does  not  readily  generalize  to  three 
dimensions. 


4.  They  require  a  priori  analytical  estimates  of  the  singularities  of  solutions  of  the  integral  equation. 

It  is  possible  to  exploit  known  analytical  estimates  characterizing  solutions  of  boundary  integral 
equations  on  singular  domains  (like  those  for  domains  with  corners  in  [7])  to  rapidly  invert  those 
equations.  However,  obtaining  even  those  results  known  today  has  been  a  laborious  process  re¬ 
quiring  great  effort,  and  the  current  state  of  affairs  is  far  from  satisfactory.  Estimates  for  planar 
domains  are  not  yet  comprehensive  and  results  for  surfaces  are  almost  entirely  lacking.  More¬ 
over,  this  approach  suffers  from  the  same  fundamental  problem  as  classical  quadrature  techniques; 
namely,  it  requires  adaptation  to  individual  cases  which  makes  it  cumbersome  to  use  in  practice. 

Unlike  previous  approaches,  our  method  does  not  rely  on  classical  quadrature  techniques  and  does 
not  require  dense  meshes  of  discretization  nodes  near  corner  points,  the  modification  of  the  underlying 
boundary  integral  formulation,  or  a  priori  analytical  estimates  of  corner  singularities.  Rather,  corner 
singularities  are  characterized  numerically  (via  orthonormal  bases)  and  efficient  quadrature  rules  which 
obviate  the  need  for  dense  meshes  of  discretization  nodes  are  precomputed  using  those  characterizations. 
The  scheme  generalizes  easily  to  a  wide  variety  of  boundary  value  problems  as  well  as  to  to  the  case  of 
two-dimensional  surfaces  with  singularities. 

Our  approach  is  most  similar  to  the  scheme  proposed  in  [9],  which  appears  to  be  the  first  paper 
describing  the  solution  of  large-scale  problems  on  domains  with  corners,  in  which  dense  meshes  of  dis¬ 
cretization  nodes  are  coupled  with  a  compression  scheme  for  the  resulting  linear  systems.  The  principal 
differences  between  the  Helsing-Ojala  algorithm  and  our  scheme  are:  (1)  our  formalism  (quadrature 
rules)  allows  for  “compression”  to  be  performed  a  priori  at  the  time  the  quadrature  rules  are  con¬ 
structed  rather  than  on-the-fly  for  each  specific  problem  as  in  [9],  and  (2)  the  Helsing-Ojala  scheme 
involves  constructs  which  do  not  appear  to  readily  generalize  to  the  case  of  surfaces  with  singularities. 

This  paper  is  divided  into  seven  sections.  After  dispensing  with  preliminary  results  pertaining  to 
quadrature  and  interpolation  in  Section  2,  we  describe  a  very  general  Nystrom  framework  for  the 
discretization  of  integral  equations  in  Section  3.  In  Section  4,  we  review  the  principal  tool  of  [4],  a 
procedure  for  the  construction  of  a  basis  spanning  the  restriction  of  solutions  of  a  particular  boundary 
integral  equation  to  a  curve  segment.  The  construction  of  quadratures  for  the  efficient  discretization 
of  a  single  given  curve  segment  is  described  in  Section  5,  while  Section  6  describes  a  procedure  for  the 
construction  of  a  collection  of  universal  quadratures  for  a  very  general  class  of  planar  domains  with 
corner  points.  Finally,  numerical  results  are  reported  in  Section  7. 


2.  Generalized  quadrature  and  interpolation 

In  this  section,  when  X  is  an  n  x  to  matrix,  we  will  denote  by  o'j(X)  its  jth  largest  singular  value. 
Moreover,  for  j  >  min (n,m)  we  define  a j(X)  =  0. 

2.1.  Discretization  of  square  integrable  functions.  We  shall  say  that  a  quadrature  rule  with  nodes 
X\, . . .  ,xn  £  [a,  6]  and  positive  weights  w i, . . . ,  wn  discretizes  a  collection  of  square  integrable  functions 
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fx,. . .  ,fm  defined  on  the  interval  [a,  b]  if 


fi(x)fj(x)dx 


^2 

1=1 


holds  for  all  z  =  1, . . . ,  to  and  j  =  1, . . . ,  to. 

If  x\, . . . ,  xn,  w  i, . . . ,  wn  is  a  quadrature  discretizing  a  collection  of  functions  in  L'2([ci,  6]), 

then  the  map  T  from  the  span  S  of  the  fj  to  the  Euclidean  space  R"  taking  the  function  /  to  the  vector 


/  f(Xl)y/Wl  \ 
f(x  2)v/W2 

\  f(xn)y/Wn  J 


(2.1) 


is  a  Hilbert  space  isomorphism  of  the  subspace  S  onto  the  subspace  of  the  Euclidean  space  R”  spanned 
by  the  vectors 


/  \ 
fi(x2)V^ 


(  fm(x \)^/Wl  \ 

fm(x  2)V/^2 


\  fl(xn)^/Wn  )  \  Sm{xn)y/Wn  ) 


2.2.  Numerical  rank.  The  numerical  rank  of  a  matrix  A  to  precision  e  is  defined  to  be  the  least 
integer  k  such  that  ak+i(A)  <  e.  Moreover,  we  define  the  numerical  rank  to  precision  e  of  a  collection 
fii . . . ,  fm  of  square  integrable  functions  on  the  interval  [a,  b]  to  be  the  numerical  rank  to  precision  e  of 
the  matrix 

/  fl(Xl)^/Wl  f2{xl)^/Wl  •••  fm(Xl)y/Wl  \ 

fl{X2)^/W2  f-2{X2)y/W2  ■■■  fm(x  2)\/W2 

•  .  •  5 

\  .fl(xn)^Wn  f2  {x^y/Wn  ...  fm(xn)y/w^  J 

where  x\, . . . ,  xn,  w i, . . . ,  wn  is  any  quadrature  discretizing  /i, . . . ,  fm. 


2.3.  Rank- revealing  QR  decompositions.  A  partial  QR  decomposition  for  an  to  x  n  matrix  A  with 
m  >  n  is  a  factorization  of  the  form 

An  =  (  Ro  £)■  (“) 

where  Q  is  an  to x to  orthogonal  matrix,  n  is  an  n x n  permutation  matrix,  f?n  is  a  kxk  upper  triangular 
matrix  with  nonnegative  diagonal  entries,  R\2  is  a  k  x  (n  —  k )  matrix,  and  I?22  is  &n  (zn  —  fc)  x  (n  —  k) 
matrix.  For  any  factorization  of  this  type,  we  have 

o’* (J?n)  <  CTj(A)  and  crj(R22)  >  ak+j(A) 

for  1  <  i  <  k  and  1  <  j  <  (n  —  k).  In  [8],  it  is  shown  that  for  any  to  x  n  matrix  A  with  m  >  n  and  any 
integer  1  <  k  <  n,  there  exists  a  factorization  of  this  form  such  that  ||J?)"11J?i2||oo  <  d(k,n), 

&i{A)  <  c(k,n)  <Ti(Ru) 


holds  for  i  =  1, . . . ,  k,  and 

c(k,n)  ak+j(A)  >  aj{R22) 

holds  for  j  =  1  k,  where  c(k,n)  =  yjl  +  k(n  —  k)  and  d(k,n)  =  1.  Moreover,  [8]  gives  a 

stable  algorithm  for  the  computation  of  partial  QR  factorizations  satisfying  bounds  of  this  type  with 
+  nk(n  —  k)  and  d  =  y/n  which  requires  0(mn2)  floating  point  operations  —  that  is,  the  same 
asymptotic  complexity  as  the  well-known  pivoted  Gram-Schmidt  algorithm2. 


2In  fact,  the  performance  of  the  scheme  described  in  [8]  depends  on  a  parameter  which  affects  both  its  running  time 
and  the  resulting  bounds.  The  values  reported  here  represent  but  one  possible  configuration. 
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Remark  2.1.  Our  implementation  of  the  algorithm  described  in  this  paper  uses  the  pivoted  Gram- 
Schmidt  procedure  with  reorthogonalization  in  place  of  the  algorithm  of  [8]  to  compute  rank-revealing 
QR  factorizations.  As  is  observed  in  [8]  and  discussed  thoroughly  in  the  monograph  [2],  the  pivoted 
Gram-Schmidt  algorithm  with  reorthogonalization,  which  is  easier  to  implement  and  somewhat  faster 
than  the  algorithm  of  [8] ,  works  well  in  practice  despite  the  existence  of  counterexamples  which  show 
that  it  can  fail  in  certain  circumstances. 


2.4.  Generalized  Chebyshev  quadratures.  A  quadrature  formula  will  be  referred  to  as  a  Chebyshev 
quadrature  for  a  set  of  2 n  linearly  independent  functions  <pi, . . . ,  02n  :  [a,  6]  — >  M  if  it  consists  of  2 n 
nodes  and  2n  weights  and  integrates  f>i,  for  all*  =  1, ... ,  2  n. 

The  construction  of  a  Chebyshev  quadrature  for  an  orthonormal  collection  of  square  integrable 
functions  u\, . . . ,  Uk  defined  on  an  interval  [a,  b }  is  trivial  given  a  preexisting  quadrature 


Xi,  .  .  .  ,  Xn ,  W\ , . . . ,  wn 

integrating  products  of  those  functions.  Let 


U  = 


/  Ul(Xl)y/wf 

Ul{x2)s/W2  •• 

■  'Ul(xn) 

ywf.  \ 

U2{x1)^/w1 

U2(X2)^/W2  ' 

■  u2{xn) 

\  uk{xi)y/W[ 

Uk(x2)^/W2  ' 

•  Uk(xn) 

s/wa  y 

and 


/  n\ 

r-2 

V  rfc  y 


where  r»,  *  =  1,, . . . ,  k,  is  defined  by 


fb 

r.i  =  /  u.i(x)dx. 
J  a 


By  virtue  of  the  orthonormality  of  the  u\ , . . . ,  Uk  and  the  requirement  that  the  quadrature  integrates 
products  of  the  iq,  the  rows  of  the  matrix  U  are  orthonormal.  It  follows  that  the  matrix  U  has  k 
nonzero  singular  values,  all  of  which  are  1.  The  results  of  [8]  discussed  in  Section  2.3  now  imply  that  a 
decomposition  of  U  of  the  form 

un  =  Q  (  R\\  R\2  ) , 

where  Q  is  an  orthogonal  k  x  k  matrix,  R\2  is  an  k  x  (n  —  k)  matrix  and  f?n  is  a  k  x  k  matrix  such 
that 

/1-L  =  n  <  <rj{Rn)  <  1  (2-3) 

i/l  +  nk(n  —  k) 

for  1  <  j  <  k,  can  be  computed  stably  in  at  most  0(n3)  operations.  A  vector  z  with  at  most  k  nonzero 
entries  such  that 


Uz  = 


(  ri\ 

r-2 

V  rfc  y 


can  be  computed  by  solving  the  k  x  k  linear  system 

Rnz  =  Q*r, 

the  condition  number  of  which  is  bounded  by 

+  nk(n  —  k ) 


(2.4) 
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by  virtue  of  the  inequality  (2.3),  and  setting  2  to  be 


ir 


l~z\ 

0 


w 


If  *i, . . .  ,ik  denote  the  indices  of  the  nonzero  components  of  z  and  we  let 


Uj  =  Xi .  and  Vj  =  z ^  ■S/Wij 

for  j  =  1, ...  ,k,  then  we  have 


/  ui(yi) 

^1(2/2)  • 

■  «i  (yk)  \ 

(  Vl  \ 

u2(yi) 

^2(2/2)  • 

•  u2{yk) 

v2 

\  uk{yi) 

uk{y2)  ■ 

■  uk(yk )  / 

\  vk  J 

that  is,  yi, . . .  ,yk,vi, . . .  ,Vk  is  a  Chebyshev  quadrature  for  u±, ...  ,uk. 


2.5.  Chebyshev  quadratures  and  interpolation.  If  f  is  a  linear  combination 


fc 

fix)  =Z«JuJ(x) 

3=1 


(2.5) 


of  the  orthonormal  functions  u-\ , . . . ,  uk  in  L2[a ,  5]  and  x\, . . . ,  xn,  «q, . . . ,  wn  is  a  quadrature  integrating 
products  of  the  the  Uj,  then  the  coefficients  oq, . . . ,  a*  in  the  expansion  (2.5)  can  be  computed  stably 
from  the  values  of  the  function  /  at  the  quadrature  nodes  Xi, ...  ,xn.  In  particular,  let  U  be  the  n  x  k 
matrix  whose  entries  are  given  by 

Uij  Uj  (x i  )  yj Wj , 


let 


F  = 


(  1)  \ 

/(*  2) 


V  f{Xn)  ) 


and  denote  by  W  the  n  x  n  diagonal  matrix  with  entries 


W»  = 


Then 


f:) 


U*WF. 


V  ak  / 


Note  that  UXJ*(WF)  =  WF  since  WF  is  in  the  span  of  the  column  space  of  U  by  assumption. 
The  computation  of  oq, . . .  ,aq.  in  this  fashion  is  entirely  numerically  stable  since  the  rows  of  U*  are 
orthonormal  and  IT  is  a  diagonal  matrix.  If  a  scheme  is  available  for  evaluating  the  functions  Uj  at 
arbitrary  points  x,  then  this  mechanism  provides  a  stable  means  for  interpolating  functions  in  the  span 
of  «i, . . . , Uk  from  Xi, ...  ,xn  to  arbitrary  points. 

In  fact,  interpolation  from  the  nodes  of  a  Chebyshev  quadrature  for  U\, . . .  ,Uk  constructed  in  the 
manner  described  in  the  preceding  section  is  also  stable.  To  see  this,  suppose  Xjt, . . . ,  xik  are  the  nodes 
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of  such  a  quadrature.  Then,  by  virtue  of  the  inequality  (2.3)  appearing  in  the  preceding  section,  the 
condition  number  of  the  matrix  U  defined  by 


/  «l(Zu) 

\/Wil 

U2(xil) 

■ 

•  Uk(Xi J 

\ 

u  = 

Ul(xi2) 

s/Wiz 

u2{xi2) 

s/m?  ■ 

•  Uk{xi2) 

\fWi2 

\  Mxik) 

U2{xik) 

■  uk(xik) 

/ 

is  bounded  by  -\/l  +  nk(n  —  k).  By  construction,  the  vector 


/(s*l) 

\ 

f(Xi2) 

/ 

is  in  the  span  of  the  columns  of  U,  so  it  follows  that  the  linear  system 


(  ai  ) 

(  fixijy/w^  \ 

a2 

f(xi2)^Ju- 

U 

— 

\  ak  ) 

\  /(^iJv7^  / 

admits  a  unique  solution  which  can  be  stably  computed. 


3.  Nystrom  discretization 


The  Nystrom  discretization  of  the  integral  equation 

Acr(a;)  +  J  K(x,y)a(y)dS(y)  =  u(x)  (3.1) 

operates  by  fixing  local  representations  for  solutions  a  of  the  form 

N 

°{y)  -^oijcfijiy) 
j= i 

and  approximating  the  integral  in  Equation  (3.1)  via  quadrature  formulae  integrating  functions  of  the 
form 

N 

K(x,y)cij  (j>j  (y) .  (3.2) 

3= 1 

Typically,  more  than  one  quadrature  formula  of  this  type  is  required,  with  each  integrating  functions 
of  the  form  (3.2)  for  various  values  or  ranges  of  values  of  x.  See  [1]  for  a  detailed  discussion  of  the 
Nystrom  discretization  of  boundary  integral  equations. 

We  now  describe  a  very  general  Nystrom  framework  for  the  discretization  of  Laplace  and  Helmholtz 
boundary  integral  equations  of  the  form  (3.1).  We  begin  by  assuming  that  T  is  divided  into  n  curve 
segments,  Tj  , . . . ,  Tn,  not  necessarily  of  equal  length.  For  each  curve  segment  Tj  we  will  require  the 
following: 

(1)  An  orthonormal  collection  of  k  basis  functions  </>i, . . .  ,<pk  in  L2(Tj), 

(2)  A  scheme  for  interpolating  the  basis  functions  from  nodes  Ai , . . . ,  Afc  to  their  values  at 

arbitrary  points, 

(3)  A  “far”  quadrature  formula  exact  for  integrals  of  the  form 


K(x,y)ai(y)  dS(y), 


i  =  1, . . . ,  k,  when  x  sufficiently  distant  from  T;, 
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(4)  A  “near”  quadrature  formula  exact  for  integrals  of  the  form 


K{x1  y)&i(y)  dS(y), 


i  =  1, . . . ,  k,  when  x  is  outside  of,  but  close  to,  F;/ , 

(5)  A  set  of  k  “diagonal”  quadrature  formulas  each  of  the  form 

r  Mi 

/  K{X,y)a{y)dS{y)  «  K(Xi,  yi)a(yi)wi, 

Jri  i=i 


exact  whenever  cr  is  in  the  span  of  cr\ , . . . ,  cr*, . 

The  method  proceeds  under  the  assumption  that  the  restriction  of  the  unknown  solution  a  in  equation 
(3.1)  to  the  curve  segment  Tj  can  be  represented  as  linear  combinations  of  the  basis  functions  <f>  i, ...  ,<f>k 
for  Tj.  Let  Fj  and  r,;  be  two  curve  segments,  not  necessarily  distinct.  We  will  denote  by  <j>i, . . . ,  <f>n  the 
basis  functions  on  the  curve  segment  F7;  by  si,..  „sn  the  interpolation  nodes  on  the  curve  segment 
Fy ;  and  by  ti, . . .  ,tm  the  interpolation  nodes  on  the  curve  segment  I',.  The  integral  equation 

Tija(x )  =  u(x),  (3.3) 

where  T)?-  is  the  integral  operator  mapping  functions  on  Tj  to  functions  on  Tj  defined  by 

Tijcr(x)  =  f  K(x,y)a(y)dy  (3.4) 

•W 

is  then  discretized  by  repeating  the  following  sequence  of  steps  for  each  interpolation  node  t  on 

(1)  The  appropriate  quadrature  formula  X\ , . . .  ,Xi,Wi, . . .  ,Wi  for  functions  of  the  form 

K(t,  s)<ju(s),  u=  1, . . .  ,n, 

is  determined.  That  is,  depending  on  the  location  of  t  relative  to  the  curve  segment  F;7 ,  either 
the  “far”  quadrature  rule,  the  “near”  quadrature  rule,  or  one  of  the  diagonal  quadrature  rules  is 
selected. 

(2)  The  kernel  K(x,y)  is  evaluated  at  the  points  (t,xr)  for  r  =  1 and  the  1  x  l  vector  v  with 
entries 

vr  =  K(t,  xr)wr 

is  formed. 

(3)  The  1  x  l  vector  v  is  multiplied  on  the  right  by  the  l  x  n  matrix  interpolating  the  basis  functions 
4>i ,<t>n  from  the  interpolation  nodes  Si, ... . . ,  sn  on  Tj  to  the  quadrature  nodes  x\ ,xi. 

(4)  The  entries  {as}  of  the  resulting  lxn  vector  give  a  single  linear  equation 

a1a(s1)  +  a2a(s2)  +  ■■  ■  ana(sn)  =  u(t) 
constraining  the  values  of  the  solution  cr  at  the  nodes  si, . . . ,  sn. 

The  result  of  repeating  this  procedure  for  each  of  the  m  interpolation  nodes  is  an  m  x  n  linear  system 
of  the  form 


an 

ai2 

n 

) 

(  o-(si)  \ 

(  u(ti) 

\ 

an 

«22 

&2  n 

cr(s2) 

— 

u(t2) 

ami 

am2 

•  ® mn 

) 

K  a(sn )  ) 

\  u(tm) 

/ 

discretizing  the  integral  equation 

TijO  =  u. 

Repeating  the  above  procedure  for  each  pair  of  curve  segments  Fj  and  Tj,  and  accounting  for  the 
constant  term  in  equation  (3.1)  results  in  a  discrete  system  of  N  equations  in  N  unknowns  of  the  form 

Xx  +  Ax  =  y, 


where  N  is  equal  to  the  sum  of  the  number  of  interpolation  nodes  Uj  on  each  curve  segment  Fj  over 
j  =  1, ...  ,n  and  A  is  a  matrix  formed  by  concatenating  the  discrete  matrices  representing  the  Tij. 


Solving  the  amalgamated  system  yields  the  values  of  the  unknown  function  a  in  (3.1)  at  the  inter¬ 
polation  nodes  of  each  of  the  curve  segments  I’r  The  value  of  a  at  any  point  x  on  T  can  then  be 
computed  in  0(1)  operations  using  the  appropriate  interpolation  formula.  Moreover,  the  value  of  a 
layer  potential 

u(x)  =  J  K(x,y)a(y)dy 

can  be  computed  for  any  x  sufficiently  far  enough  away  from  the  curve  T  in  O(n)  operation  using  the 
far  quadrature  formulas  for  the  curve  segments  Tj,  j  =  1  For  points  close  to  the  curve,  an 

adaptive  Gaussian  quadrature  scheme  which  relies  on  the  ability  to  evaluate  the  charge  distribution  at 
any  point  via  interpolation  can  be  used  to  compute  the  value  of  the  layer  potential. 

4.  Charge  bases 

In  this  section,  we  show  that  under  very  mild  assumptions,  a  small  finite  orthonormal  basis  spanning 
the  restrictions  of  solutions  of  a  Laplace  or  Helmholtz  boundary  integral  equation  on  a  contour  T  to  a 
small  curve  segment  r0  C  T  can  be  constructed.  We  will  refer  to  such  a  basis  as  a  “charge  basis”  for 
the  contour  To  (the  underlying  integral  equation  will  always  be  apparent  from  context). 

In  the  interests  of  brevity,  we  will  restrict  our  attention  here  to  the  boundary  integral  equation 

-\a(x)  +  Jr  {)±v(y)  ■  vy  +  H°  \x  -  y\)  cr(y)  dS(y) =  u(x )>  (4-1) 

where  Hq  is  the  Hankel  function  of  zeroth  order  and  v{y)  denotes  the  outward  unit  normal  vector  on  the 
contour  T.  The  equation  (4.1)  is  typically  used  to  obtain  a  solution  of  the  exterior  Dirichlet  boundary 
value  problem  for  the  Helmholtz  equation  with  wavenumber  k  (see,  for  instance,  [5]).  The  argument 
given  here  clearly  applies  to  a  variety  of  other  boundary  integral  equations  with  minor  modifications. 

4.1.  Partial  wave  expansions.  If  the  function  <f>  :  R2  — >  R.  satisfies  the  Helmholtz  equation 

(V2  +  k2)  (j)  =  0 

in  a  disc  D  of  radius  R  >  0  centered  at  xq  and  the  radiation  condition 

lim  <j>(tx)  e"ifct'N  =o(t~1/2) 

t—>  OO  V  / 

uniformly  at  infinity,  then  <j>  can  be  represented  uniquely  as  a  J-expansion  in  D\  i.e.,  in  the  form 

OO 

4>{x)  =  22  amJm(kr)eim0 ,  (4.2) 

m=— oo 

where  (r,  6)  denote  the  usual  polar  coordinates  centered  at  Xq  and  Jm  is  the  Bessel  function  of  the  first 
kind  of  order  m.  The  following  theorem,  which  appears  in  [10],  establishes  the  convergence  rate  of  the 
expansion  (4.2). 

THEOREM  4.1.  If  Di  C  D  is  a  disc  of  radius  R\  <  R  centered  at  Xq,  then  there  exists  a  constant 
c  >  0  such  that  for  x  £  D i  and  N  >  |fc|  R\, 

N 

H X )  -  22  amJm{kr)em6 

m——N 

4.2.  Rank  of  interaction.  If  To  and  Ti  are  two  disjoint  compact  connected  Lipschitz  curves  in  the 
plane  and  K{x,  y)  is  one  of  the  potential  theoretic  kernels  considered  in  this  paper,  then  the  mapping 
defined  by 

Tf{x)=  f  K(x,y)f(y)dS(y)  (4.3) 

■M 

is  compact  as  an  operator  L 2  (To)  — >  L2  (Ti).  We  shall  define  the  rank  of  interaction  of  the  curve 
segments  T 0  and  T  i  under  the  kernel  K  as  the  least  integer  m  such  that  either  m  —  1  singular  values 
of  T  are  less  than  the  prescribed  precision  e. 
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Figure  1:  The  rank  of  interaction  of  the  two  curve  segments  To  and  Ti  under  the  potential  theoretic 
kernels  considered  in  this  paper  is  low  regardless  of  the  angle  of  the  corner. 


If  the  rank  of  interaction  of  the  operator  (4.3)  is  m  to  precision  e,  then  there  exists  a  quadrature 
formula  of  the  form 


K(x,  y)f(y)dS(y)  «  ^ K(x,yJ)f{yj)w:j , 


with  the  yj  lying  in  Ti,  which  holds  with  precision  approximately  e  for  a;  £  To  (note  that  the  kernels 
considered  in  this  paper  are  bounded  almost  everywhere  for  leTj  and  y  €  Id).  A  discussion  of  the 
existence  and  numerical  computation  of  such  quadratures  can  be  found  in  [4]. 

Although  the  rank  of  interaction  between  two  curve  segments  under  one  of  the  potential-theoretic 
kernels  considered  here  can  be  high,  for  most  curves  appearing  in  practice  the  singular  values  of  the 
integral  operators  corresponding  to  interactions  between  neighboring  curve  segments  decay  exponen¬ 
tially  (see,  for  instance,  [6],  which  contains  several  relevant  examples).  In  the  particular  case  of  corner 
regions,  which  are  the  focus  of  this  paper,  the  rank  of  interaction  between  a  neighborhood  of  a  corner 
and  the  neighboring  portions  of  the  curve  is  typically  low  regardless  of  the  angle  of  the  corner.  The 
situation  is  illustrated  in  Figure  1. 


4.3.  Construction  of  charge  bases.  In  what  follows,  we  shall  fix  a  simply-connected  domain  O 
in  the  plane  whose  boundary  T  is  a  compact,  connected  Lipschitz  curve.  Our  aim  is  to  produce  an 
orthonormal  basis  approximately  spanning  the  set  of  solutions  of  the  boundary  integral  equation  (4.1) 
to  a  small  curve  segment  Tq  C  T.  We  will  denote  by  K(x,y )  the  integral  kernel 


(ylv(v)  '  vy  +  ^  Ho  (k  \x  -  y\) 

appearing  in  equation  (4.1).  Moreover,  will  shall  let  B  be  a  disc  of  minimum  radius  which  contains 
To  and  we  will  denote  by  Tj  the  portion  of  the  contour  T  contained  in  2 B  \  B.  Finally,  we  will  let  T2 
denote  the  portion  of  the  curve  T  contained  in  the  complement  of  the  ball  2 B.  Figure  2  depicts  the 
situation. 

The  boundary  integral  equation  (4.1)  can  be  rearranged  as 


i(T(x)+  I  K(x,y)cr(y)dS{y)  =  u{x)  - 
z  Jr„ 


/  K(x,y)cr(y)dS(y)  —  /  K(x,y)a(y)dS{y). 
/  r2  Jr1 


Under  the  assumption  that  the  right  hand  side  u(x)  satisfies  the  Helmholtz  equation  in  the  disc  2B, 
we  can,  by  virtue  of  Theorem  4.1,  introduce  the  approximation 


u{x) 


+  [  K(x,y)a(y)dS(y) 
Jr2 


N 

Y,  amJm{kr)eimd , 

m=—N 


(4.5) 
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Figure  2 


which  holds  for  x  £  To-  Moreover,  as  described  in  Section  4.2,  when  a:  €  To  the  third  term  on  the  right 
hand  side  of  equation  (4.4)  can  be  approximated  as 


.  M 

/  K(x,y)f(y)d,S(y)  «  ^ ^K(x,yj)f(yj)wj , 

JT'  3=1 

where  the  yj  lie  in  F  |  and  M  is  the  rank  of  interaction  of  To  with  Ti .  It  follows  that  the  restriction  of 
< t  to  the  curve  segment  To  satisfies  the  integral  equation 

N  M 


1 

~2a 


(x)+f  K(x,y)a(y)dS(y)  =  ^  amJm(kr)eime  +  y^iK(x,yj)f(yj)wj.  (4.6) 

Jr0  rn—  —  N  j=  1 


It  is  now  clear  how  to  form  a  charge  basis  for  To-  We  observe  that  the  restriction  of  (4.1)  to  To  is 
invertible  as  an  operator  L2{Tq)  — >  L2(To)  and  form  the  collection  of  functions  obtained  by  solving  the 
restricted  integral  equation  for  each  of  the  functions  of  the  form 

Jm(kr)elrne  and  K(x,yj) 

appearing  in  (4.6).  The  resulting  functions  are  orthonormalized  in  order  to  form  a  charge  basis. 


4.4.  Convergence  estimate.  In  this  section,  we  derive  a  crude  estimate  on  the  number  of  charge 
basis  functions  required  to  achieve  a  given  L°°(ro)  precision  e  for  the  approximation  (4.6).  And  of 
course,  this  quantity  is  related  to  the  L2(To)  error  in  the  approximation  of  restricted  solutions  by  the 
charge  basis  through  the  L2(To)  norm  of  the  inverse  of  the  restriction  of  the  integral  operator  (4.1)  and 
the  arclength  of  the  curve  To- 

In  particular,  if  N  is  chosen  so  that 


N  >  max 


-  log  (e)  +  log  (c)  \ 
log  (2)  )' 


(4.7) 


where  R\  is  the  radius  of  B  and  c  is  the  constant  appearing  in  Theorem  4.1,  then  the  approximation  (4.5) 
is  guaranteed  to  achieve  precision  e.  The  rank  of  interaction  M  of  To  and  Ti  depends  on  the  geometry 
of  the  contour  T  and  the  wavenumber  k ;  however,  in  typical  cases  —  especially  when  T0  is  taken  to  be 
a  small  region  around  a  corner  point  —  M  behaves  in  a  similar  fashion. 

Thus  the  number  of  charge  basis  functions  required  for  precision  e  is  roughly 


max  (|fc|  i?i,  0(—  log  (e))). 


In  other  words,  the  convergence  of  approximations  with  respect  to  charge  bases  is  exponential  once  a 
certain  number  of  “terms  per  wavelength”  has  been  achieved.  This  means  that  the  number  of  basis 
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functions  necessary  to  represent  a  charge  distribution  on  many  types  of  complicated  curve  segments 
(e.g.  corner  regions)  is  comparable  with  the  number  required  to  represent  a  charge  distribution  on  a 
smooth  curve  segment  using  typical  approaches  (e.g.,  piecewise  Gaussian  polynomials).  This  estimate 
is  consistent  with  our  numerical  experiments,  which  are  presented  in  Section  7. 

Remark  4.1.  This  estimate  only  provides  a  very  crude  upper  bound  for  the  dimension  of  the  resulting 
charge  basis;  charge  bases  are  almost  always  much  smaller  than  this  estimate  suggests.  Moreover,  in 
practice,  the  number  of  charge  basis  functions  seems  to  have  only  a  weak  dependence  on  the  rank  of 
interaction  of  To  andT i. 

Remark  4.2.  In  the  case  of  the  boundary  integrals  related  to  Laplace’s  equation,  considerations  in¬ 
volving  the  wavenumber  are  (obviously)  unnecessary.  For  those  integral  equations,  in  typical  cases  the 
number  of  charge  basis  functions  behaves  as  0(—  log(e)). 

5.  Purpose-made  quadratures 

In  this  section,  we  describe  an  algorithm  for  the  construction  of  a  charge  basis  cr\ , . . . ,  cr^  for  a 
Laplace  or  Helmholtz  boundary  integral  equation 

\cr(x)  +  J  K(x,y)a(y)  dS(y)  =  u(x)  (5.1) 

over  a  corner  curve  segment  To  C  T  as  well  as  an  algorithm  for  the  computation  of  purpose-made 
quadratures  for  its  discretization  over  r0  given  a  charge  basis  for  T0. 

5.1.  Construction  of  a  charge  basis.  The  algorithm  of  this  section  takes  as  input  a  precision  e,  a 
parameterization  r  :  [—  1, 1]  — >  T0  that  maps  0  to  the  corner  point  of  T0,  and  two  integer  parameters  s 
and  k  whose  roles  will  be  described  below.  It  proceeds  as  follows: 

1.  Form  a  simply-graded  mesh1  on  the  corner  curve  segment  To  by  constructing  the  subintervals  of 
[—1, 1]  with  endpoints 

—  and - -  for  j  =  1,  2, . . . ,  s, 

V  V  J 

and  taking  the  image  of  the  resulting  subintervals  under  the  parameterization  r. 

2.  Discretize  the  restricted  boundary  integral  equation 

Acr(x)  +  f  K(x,y)a(y)  dS(y)  =  u(x)  {x  €  r0)  (5.2) 

Jt0 

via  the  Nystrom  method  under  the  assumption  that  the  solution  a  can  be  represented  as  a  kth 
order  piecewise  Legendre  expansion  over  the  simply-graded  mesh  formed  in  Step  1.  Denote  by 
Xi, . . . ,  xm,  w  i, . . . ,  wm  the  nodes  and  weights  of  the  piecewise  Legendre  quadrature  which  represents 
solutions  of  (5.2)  on  T0. 

3.  Solve  the  discrete  linear  system  formed  in  Step  2  for  a  collection  of  right-hand-sides  consisting  of 
multipoles  and  functions  of  the  form 

K(x,y3), 

where  the  y3  are  points  on  T  adjacent  to  T0.  The  order  of  the  multipole  functions  is  chosen  so 
as  to  ensure  very  high  accuracy  in  the  approximation  of  potential  functions  on  To-  Similarly,  the 
points  y:j  adjacent  to  To  are  chosen  to  be  the  nodes  of  a  quadrature  which  discretizes  the  operator 
T  :  L2(r0)  ->  L2(T i)  given  by 

Tf( x)=  f  K(x,y)f{y)dy 
Jr0 

to  high  accuracy. 

Denote  the  resulting  functions  by  Ui, . . .  ,un.  Note  that  the  values  of  the  Uj  at  arbitrary  points  on 
r0  can  be  calculated  using  standard  interpolation  techniques  for  piecewise  Legendre  polynomials. 

1Tho  terminology  is  adapted  from  [9].  Also  see  [4]. 
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4.  Form  the  in  x  n  matrix  A  whose  entries  are  given  by 

Aij  =  Uj(xi)^w~i. 

5.  Form  a  rank-revealing  QR  decomposition 

'4n  =  0(*1  £() 

for  A  and  let  k  be  the  least  integer  such  that  the  fcth  singular  value  of  A  is  less  than  the  input 
precision  e. 

6.  Define  the  charge  basis  functions  a i, . . . ,  oy,  via  the  formula 

(7 j  (Xj  )  —  Qij  /  sJ'Wi , 

where  Qij  refers  to  the  ij th  entry  of  the  matrix  Q  formed  in  the  preceding  step.  Note  that  since 
the  values  of  the  a j  are  known  at  the  points  of  a  piecewise  Legendre  quadrature,  their  values  at 
any  point  on  To  can  be  readily  obtained  via  interpolation. 

The  parameters  s  and  k  should  be  chosen  so  as  to  ensure  that  the  piecewise  Legendre  quadrature 
used  to  discretize  solutions  of  the  restricted  boundary  integral  equation  (5.2)  integrates  products  of  the 
solutions  obtained  in  Step  3. 

Remark  5.1.  Simply-graded  meshes  are  not  necessarily  sufficient  to  represent  solutions  of  the  bound¬ 
ary  integral  equations  to  double  precision,  at  least  if  the  computations  are  performed  using  double  preci¬ 
sion  arithmetic.  However,  the  precomputations  described  in  this  section  may  be  performed  in  extended 
precision  (FORTRAN  REAL*  16)  arithmetic  if  necessary  to  achieve  higher  accuracy. 

5.2.  Construction  of  the  quadrature  formulae.  In  this  section,  we  detail  the  construction  of  the 
quadrature  formulae  for  the  discretization  of  the  boundary  integral  equation  (5.1)  over  To- 

The  algorithm  of  this  section  takes  as  input  the  quadrature  X\, . . . ,  xm,  w i, . . . ,  wm  and  the  charge 
basis  <7i <Tfc  formed  using  the  procedure  of  the  preceding  section  as  well  as  an  integer  parameter 
fciege-  We  may,  by  virtue  of  the  parameterization  r  :  [—1, 1]  — ■>  To  regard  the  functions  a3  as  given  over 
the  interval  [—1,1], 

The  algorithm  proceeds  as  follows: 

1.  Construct  a  Chebyshev  quadrature  for  the  charge  basis  functions  o\, . . . ,  ay,.  The  nodes  Ai, . . . ,  A*, 
of  this  quadrature  serve  as  discretization  nodes  for  Nystrom  procedure  of  Section  3.  Note  that  the 
nodes  Ay, .. . ,  A  a,  are  guaranteed  to  be  stable  interpolation  nodes  for  the  charge  basis  functions  (as 
described  in  Section  2.5). 

2.  Form  the  “far”  quadrature  required  by  the  Nystrom  procedure  by  constructing  a  Chebyshev  quad¬ 
rature  for  functions  of  the  form 

Gj{x)  (p(x)  +q(x)) , 

where  1  <  j  <  k,  p(x )  is  a  polynomial  of  degree  k iege  on  the  interval  [—1, 0],  and  q(x )  is  a  polynomial 
of  degree  klege  on  the  interval  [0, 1]. 

3.  Construct  the  “near”  quadrature  by  constructing  a  Chebyshev  quadrature  for  functions  of  the  form 

log \t  -  x\  C Tj(x), 

where  j  =  1, ...  ,k  and  t  is  a  point  on  T  close  to  To- 

4.  For  each  i  =  1, . . . ,  k,  construct  a  diagonal  quadrature  for  functions  of  the  form 

K(\i,x)<Tj(x), 

j  =  1,  •  ■  • ,  k. 

5.  For  each  of  the  k  +  2  quadratures  constructed  in  the  last  section,  form  the  matrix  interpolating  the 
charge  basis  functions  from  the  discretization  nodes  Ai, . . . ,  A*,  to  the  quadrature  nodes.  This  can 
be  achieved  via  the  procedure  described  in  Section  2.5. 
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6.  Universal  quadratures  for  domains  with  corners 


We  now  indicate  how  the  procedure  of  the  preceding  section  can  be  modified  in  order  to  produce  a 
collection  of  quadrature  formulae  suitable  for  the  discretization  of  a  very  general  class  of  domains  with 
corner  points. 

We  shall  say  that  a  contour  T  is  a  corner  region  with  corner  point  at  ( Xo,yo )  if  it  admits  a  parame¬ 
terization  r(t)  =  (x(t),y(t))  such  that 


x(t) 

y(t) 


Xq  +  ait  +  a2t"  “I-  a3 ■ 
Xq  +  b±t  +  b2t2  +  63 +  . . . 


for  a  <  t  <  0 
for  0  <  t  <  b 


yo  T  Cit  C2t2  -t-  C3 G  T  . . . 
1/0  +  dit  +  d2t2  +  t'^  +  . . . 


for  a  <  t  <  0 
for  0  <  t  <  6, 


(6.1) 


where  the  series  converge  uniformly  for  t  ^  0,  and  (ai,Ci)  ^  (bi,di).  We  can  assume  without  loss  of 
generality  that  (6.1)  is  an  arclength  parameterization;  i.e.,  that  x'(t)2  +  y'(t)2  =  1.  For  a  contour  of 
this  type,  we  define  the  oriented  angle  6  of  the  corner  as  the  angle  of  the  counter-clockwise  rotation 
required  to  transform  the  vector 


lim 

t— >  0- 


x'(t) 

y'(t) 


a  1 

Cl 


into  the  vector 


lim 
t— +0+ 


x'{t) 
y'(t ) 


h 

d\ 


Furthermore,  we  define  the  left  and  right  instantaneous  curvatures  by 

=  tlim_  [x'(t)y"(t)  -  y'(t)x"(t)} 

=  2(aiC2  —  Cia2) 

K+  =  tlh?+  [x'^y"^  ~  2/,(*)a;,,W] 

=  2(bi d2  -  dib2). 


Given  the  angle  9  and  the  two  instantaneous  curvatures  K-  and  the  coefficients  a±,  a2,  bi,  b2 , 
Ci,  C2,  d\,  d2  in  the  parameterization  (6.1)  can  be  recovered  up  to  rotation  and  translation  —  that 
is,  the  coefficients  corresponding  to  a  contour  which  is  the  image  of  T  under  translation  and  rotation 
can  be  readily  derived  from  9 ,  k_,  and  k+.  Note  that  this  depends  on  our  assumption  that  the 
parameterization  is  with  respect  to  arclength  (which  reduces  the  number  of  free  variables  to  4). 

It  is  an  obvious  consequence  of  the  rotation  and  translation  invariance  of  Laplace  and  Helmholtz 
problems  that  a  charge  basis  for  corner  regions  whose  oriented  angles  and  instantaneous  curvatures  lie 
in  given  ranges 

9G[9U92] 

K-,K+  G  [«i,  K2\ 

and  such  that  the  terms  of  order  greater  than  2  in  the  parameterization  (6.1)  are  0  can  be  formed  by 
sampling  a  large  number  of  angles  9  and  curvatures  and  k+  in  the  desired  intervals  and  combining 
the  charge  bases  corresponding  to  corner  regions  with  each  possible  combination  of  these  parameters. 
However,  it  is  also  possible  to  remove  the  restriction  that  the  higher  order  terms  in  (6.1)  be  zero. 
To  do  so,  we  observe  that  a  perturbation  in  the  higher  order  terms  in  (6.1)  results  in  a  piecewise 
smooth  perturbation  of  Laplace  or  Helmholtz  integral  kernels  defined  over  the  contour  T.  By  artificially 
augmenting  the  basis  obtained  by  combining  charge  bases  for  various  values  of  9,  K- ,  with  functions 
of  the  form 

p{x)  +  q(x) 

where  p(x)  is  a  polynomial  over  a  given  degree  on  the  interval  [a,  0)  and  q{x)  is  a  polynomial  over  a 
given  degree  on  the  interval  (0,  6],  we  obtain  a  basis  which  approximately  spans  the  space  of  restrictions 
of  a  Laplace  or  Helmholtz  boundary  integral  equation  to  corner  contours  of  the  form  (6.1)  for  a  wide 
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class  of  corners  with  oriented  angles  9  £  [6\ ,  6*2]  and  instantaneous  curvatures  £  [k\  ,  •  In  a 

similar  manner,  a  collection  of  quadratures  for  the  discretization  of  boundary  integral  equations  over 
such  regions  —  that  is,  the  far,  near,  and  diagonal  regimes  described  in  the  preceding  section  —  can  be 
obtained  by  artificially  augmenting  the  collections  of  functions  which  serve  as  input  to  the  Chebyshev 
quadrature  procedure  with  bases  for  piecewise  smooth  functions. 


7.  Numerical  results 

The  algorithm  of  this  paper  for  the  construction  of  universal  quadratures  for  domains  with  corners 
was  implemented  in  Fortran  77  as  was  (the  obvious)  modification  of  that  algorithm  for  the  construction 
of  purpose-made  quadratures  for  specific  curve  segments.  The  resulting  code  was  complied  with  the 
Lahey/Fujitsu  Linux64  Fortran  Compiler  Release  8.10a.  Timings  were  performed  on  a  PC  with  an  Intel 
Core  i7  2.67  GHz  processor  and  12GB  of  memory  and  refer  to  wall  clock  time.  No  attempt  was  made 
to  parallelize  any  of  the  code. 

All  boundary  integral  equations  were  inverted  using  a  simple  direct  solver  which  is  O  (n2)  in  the 
number  of  discretization  nodes  n.  The  precise  algorithm  will  be  reported  at  a  later  date. 


7.1.  Purpose-made  quadratures. 


7.1.1.  A  Laplace  Neumann  problem.  In  this  numerical  experiment,  we  solved  the  exterior  Neumann 
problem 


Au(x)  =  0  for  x  £  fig 

lim  ^(x)  =  f(p)  for  p  £  <9H3,  C7’1) 

where  fl3  is  the  domain  shown  in  Figure  5c,  ^  denotes  differentiation  with  respect  to  the  outward 
pointing  normal  vector,  and  the  boundary  data  f(p)  was  taken  to  be  the  outward  normal  derivative 
on  <912  3  of  a  potential  function  u  arising  from  5  charges  placed  randomly  in  the  interior  of  123.  The 
solution  u(x)  of  (7.1)  was  represented  in  the  form 

u(x)  =  Jr:./  log  \x  -  y\ a{y)  dS(y).  (7.2) 

2t  Jpn3 

A  total  of  871  quadratures  nodes  were  used  to  discretize  the  resulting  boundary  integral  equation 

la(x)+^~[  Vo,  log  \x  —  y\  ■  v{x)a{y)  dS(y)  =  f(x).  (7.3) 

^  * n  JdQ3 

The  12  corner  regions  were  discretized  using  purpose-made  quadratures,  which  ranged  in  size  from 
17  to  21  nodes,  while  the  smooth  portions  of  the  curve  were  discretized  using  630  piecewise  Legendre 
quadrature  nodes.  Figure  3  shows  four  charge  basis  functions  obtained  in  the  process  of  constructing 
a  purpose-made  quadrature  for  one  of  the  corner  regions. 

In  order  to  access  the  accuracy  of  the  discretization,  the  difference  between  the  true  potential  function 
v  and  the  single  layer  potential 

u(x)  =  J  :  [  I°g  \x  -  y\ cr(y)  dS(y), 

Jon3 

where  cr  is  the  computed  inverse  of  equation  (7.3),  was  measured  at  a  collection  of  100  randomly  chosen 
points  in  the  exterior  of  the  domain  Q3  as  well  as  at  300  points  on  a  circle  T  enclosing  the  domain 
03.  The  largest  absolute  error  was  found  to  be  1.78  x  10-13  and  the  estimated  relative  L2(T)  error 
Ilu  —  UH2/IHI2  was  6.24  x  10-13. 
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(b) 


(c) 


(d) 


Figure  3:  Four  charge  basis  functions  used  in  the  discretization  of  a  Laplace  Neumann  boundary  integral 
equation  on  the  contour  d^.3. 


7.1.2.  A  Helmholtz  Dirichlet  problem.  In  this  example,  we  solved  the  Helmholtz  exterior  Dirichlet 
problem 

A u(x)  +  k2u(x )  =  0  for  x  € 

lim  u(x)  =  f(p)  for  x  €  fic,  (^-4 ) 

where  fi  is  the  “shark’s  fin”  domain  shown  in  Figure  4,  the  wavenumber  k  is  20,  and  f(p)  is  taken  to  be 
the  restriction  to  the  boundary  <912  of  a  acoustic  potential  u(x)  generated  by  5  random  charges  placed 
in  the  interior  of  the  domain  f2.  The  domain  fi  is  approximately  80  wavelengths  by  80  wavelengths  in 
size. 

The  solution  u(x)  of  the  problem  (7.4)  was  represented  in  the  form 

u(x)  =  Qz'(y)  •  V.y  +  1^  H0  ( k  \x  -  y\)  a(y)  dS(y),  (7.5) 

where  v{y)  denotes  the  outward  pointing  normal  vector  to  the  contour  dfl  at  the  point  y ,  Vy  denotes  the 
gradient  in  the  y  variable,  and  Hq  is  the  Hankel  function  of  the  first  kind  of  order  0.  This  representation 
leads  to  the  boundary  integral  equation 

-\aix)  +  Jgn  (^v(y) '  vy  +  ^  H°  ( k  \x  -  2/1)  cr(2 /)  dS(y)  =  /( x).  (7.6) 

This  integral  equation  was  discretized  using  a  combination  purpose-made  quadratures,  one  for  each  of 
the  3  corner  regions  of  <912,  and  30-point  piecewise  Legendre  quadratures  for  the  smooth  portions  of  the 
contour.  A  total  of  1163  quadratures  nodes  were  used;  the  purpose-made  quadratures  were  70,  68,  and 
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65  points.  The  error  in  the  representation  (7.5)  was  measured  at  100  random  points  in  the  exterior  of 


Figure  4:  The  “shark’s  fin”  domain  U  under  consideration  in  Section  7.1.2. 

the  domain  Q  as  well  as  at  300  points  on  a  circle  T  enclosing  O.  The  largest  absolute  error  was  found 
to  be  3.56  x  10-13  while  the  estimated  relative  L2(T)  error  was  found  to  be  7.67  x  10-12. 

7.2.  Universal  quadratures.  Using  the  procedures  described  in  Sections  5  and  6,  universal  quadra¬ 
tures  were  constructed  for  the  discretization  over  domains  with  corners  of  the  boundary  integral  equa¬ 
tions 

-^cr(^)  +  J-/  (Vy  log  \x  —  y\-  v{y)  +  1)  a(y)  dS(y)  =  u(x),  (7.7) 

^  ^  JdQ 

^a(x)  +  ^-[  V^logla;  -  y\  ■  v(y)a(y)  dS(y)  =  u(x),  (7.8) 

^  JdQ 

and 

-^(a:)  +  (^v(y)  •  Vy  +  1^  Ho  (k  \x  -  y\)  a(y)  dS(y )  =  u(x),  (7.9) 

where  Hq  is  the  Hankel  function  of  zeroth  order,  v{y)  denotes  the  outward  pointing  unit  normal  vector 
of  d Cl  at  the  point  y.  Several  collections  of  quadratures  were  constructed  for  the  boundary  integral 
equation  (7.9),  one  for  each  of  several  different  ranges  of  the  wavenumbers  k  in  the  interval  [.1,  20].  The 
integral  equation  (7.7)  arises  from  the  solution  of  the  exterior  Dirichlet  problem  for  Laplace’s  equation, 
(7.8)  arises  from  the  solution  of  the  exterior  Neumann  problem  for  Laplace’s  equation,  and  (7.9)  is 
associated  with  the  exterior  Dirichlet  problem  for  the  Helmholtz  equation  with  wavenumber  k  (see  [5], 
for  instance). 


^2 

Problem 

k 

N 

T 

Ep0t 

N 

T 

Epot 

Laplace  Dirichlet 

_ 

258 

0.02 

9.90xl0-12 

1320 

0.31 

4.86xl0-13 

Laplace  Neumann 

- 

351 

0.02 

3.80xl0“14 

1376 

1.45 

4.72xl0-12 

Helmholtz  Dirichlet 

0.1 

364 

0.71 

l.llxlO-13 

1512 

4.82 

1.72xl0-12 

0.5 

370 

0.74 

6.90xl0“14 

1528 

5.12 

3.55xl0“12 

1 

375 

0.79 

7.75xl0“14 

1532 

5.19 

7.87xl0-12 

2 

380 

0.83 

1.38x  10-13 

2504 

10.76 

6.24xl0-14 

5 

570 

1.35 

5.57xl0“12 

2532 

11.52 

2.33xl0-13 

10 

572 

1.49 

7.98xl0“14 

2552 

14.79 

1.75xl0-11 

20 

939 

2.59 

7.81xl0-14 

4492 

31.23 

6.75xl0“12 

Table  1:  Computational  results  obtained  for  the  Pacman  and  inkblot  domains  Hi  and  fl-2- 
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(a)  The  “Pacman”  domain  f2i. 


(b)  The  “inkblot”  domain  V-2  ■ 


(c)  The  “tank”  domain  Q3. 

Figure  5:  The  domains  under  consideration  in  Section  7.2.  The  Pacman  domain  is  approximately  4 
units  by  4  units  in  size,  while  the  dimensions  of  the  inkblot  domain  are  approximately  6  by  6.  The 
tank  domain  is  approximately  6  units  in  length  by  3  units  in  height. 


H3 

Problem 

k 

N 

T 

Epot 

Laplace  Dirichlet 

- 

1197 

0.76 

9.10xl0~12 

Laplace  Neumann 

- 

1385 

1.50 

4.62xl0^12 

Helmholtz  Dirichlet 

0.1 

1359 

8.73 

2.54xl0^13 

0.5 

1385 

9.29 

4.82xl0~13 

1 

1403 

9.61 

l.llxlO^12 

2 

1992 

14.74 

1.30xl0^13 

5 

2746 

25.49 

1.13xl0~13 

10 

2777 

26.71 

7.59xl0^13 

20 

3221 

30.94 

3.70xl0^13 

Table  2:  Computational  results  obtained  for  the  domain 


For  each  of  the  domains  shown  in  Figure  5,  the  exterior  Dirichlet  and  Neumann  problems  for  Laplace’s 
equation  and  the  exterior  Dirichlet  problem  for  the  Helmholtz  equation  at  various  wavenumbers  were 
solved.  In  each  case,  the  boundary  data  was  taken  to  be  an  electromagnetic  or  acoustic  potential 
generated  by  a  collection  of  5  random  charges  in  the  interior  of  the  domain  and  the  boundary  integral 
equation  was  discretized  using  universal  quadratures  for  corner  regions  and  piecewise  Legendre  quadra¬ 
tures  on  the  smooth  portions  of  the  curve.  The  resulting  solution  was  tested  by  comparing  the  value  of 
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the  true  potential  function  to  the  computed  solution  at  100  randomly  chosen  points  in  the  exterior  of 
the  domain.  Tables  1  and  2  present  the  results.  The  following  quantities  are  reported  for  each  problem: 

•  k  the  wavenumber  for  the  problem  (if  applicable); 

•  N  is  the  number  of  discretization  nodes  required  for  the  problem; 

•  T  is  the  wall  clock  time  required  to  solve  the  integral  equation;  and 

•  Epot  the  largest  absolute  error  measured  in  the  computed  potential. 


Figure  6:  The  domain  Q4,  formed  from  the  graph  of  two  Lipschitz  functions. 

7.3.  A  final  example.  As  a  final  example,  we  used  universal  quadratures  to  solve  an  exterior  Dirichlet 
problem  for  the  Helmholtz  equation  at  wavenumber  k  =  1  on  the  domain  Q4  shown  in  Figure  6.  The 
boundary  of  1^4  is  formed  from  two  Lipschitz  functions  and  has  40  corner  points.  The  domain  H4  is 
approximately  3  wavelengths  in  width  and  3  wavelengths  in  height.  The  boundary  data  for  the  problem 
was  taken  to  be  an  acoustic  potential  generated  by  5  charges  randomly  placed  in  the  interior  of  the 
domain. 

The  discretization  of  the  boundary  integral  equation  (7.9),  which  was  used  to  solve  the  boundary 
value  problem,  required  7794  discretization  nodes.  Between  54  and  70  discretization  nodes  were  required 
per  corner.  Inversion  of  the  resulting  discrete  linear  system  required  71.83  seconds  and  resulting  in  an 
approximate  charge  distributions  which,  when  compared  with  the  true  potential  function  at  a  collection 
of  100  randomly  placed  points  in  the  exterior  of  the  domain,  yielded  a  maximum  absolute  error  of 
2.13  x  1CT12. 

8.  Conclusions 

We  have  introduced  an  algorithm  for  the  construction  of  “universal  quadratures”  for  the  efficient 
discretization  of  Laplace  and  Helmholtz  boundary  integral  equations  over  a  very  general  class  of  domains 
with  corner  points.  Once  they  have  been  constructed,  these  quadratures  reduce  the  complexity  of  solving 
a  boundary  integral  equation  on  such  domains  to  approximately  that  of  solving  the  same  equation  over 
well-behaved  smooth  domains.  While  the  implementation  presented  here  is  for  certain  Laplace  and 
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Helmholtz  boundary  value  problems  over  planar  domains  with  corners,  the  extension  of  our  approach 
to  other  boundary  value  problems,  other  pathological  domains,  and  to  surfaces  with  singularities  is 
straightforward  and  will  be  reported  at  a  later  date. 
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